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Abstract 

Pattern formation in reaction-diffusion systems is of great importance in surface micro-patterning [Grzy- 

bowski et al. Soft Matter. 1, 1 14 (2005)], self-organization of cellular micro-organisms [Schulz et al. 

'f-f-> ' . Annu. Rev. Microbiol. 55, 105 (2001)] and in developmental biology [Barkai et al. FEES Journal 276, 

c/3 \ 1196 (2009)]. In this work, we apply the Lattice Boltzmann method (LBM) to study pattern formation in 

Cd ■ reaction-diffusion systems. As a first methodological step, we consider the case of a single species un- 

B' 

I ■ dergoing transformation reaction and diffusion. In this case, we perform a third-order Chapman-Enskog 

fi \ multiscale expansion and study the dependence of the Lattice Boltzmann truncation eiTor on the diffusion 

' coefficient and the reaction rate. These findings are in good agreement with numerical simulations. Fur- 

thermore, taking the Gray-Scott model as a prominent example, we provide evidence for the maturity of the 
0^ , LBM in studying pattern formation in non-linear reaction-diffusion systems. For this purpose, we perform 



o 



linear stability analysis of the Gray-Scott model and determine the relevant parameter range for pattern 



t^^ ■ formation. Lattice Boltzmann simulations allow not only to test the validity of the linear stability phase 

o ; 

O . diagram including Turing and Hopf instabilities, but also permit going beyond the linear stability regime, 

where large perturbations give rise to interesting dynamical behavior such as the so called self replicating 

/\ ' spots. We also show that the length scale of the patterns may be tuned by rescaling all relevant diffusion 

coefficients in the system with the same factor while letting all the reaction constants unchanged. 



I. INTRODUCTION 



Spatially and/or temporally var ying patterns have been observed in a variety of physical yj, |21] , 



chemical ||3|-l5|] and biological II6I-IIIII systems operating far from equilibrium. The interest in 
understanding the physics of pattern formation in these systems has been increasing steadily over 
the last few years especially after the experimental verification of Turing's idea lll2ll . In chemical 
and biological systems for instance, macroscopic reaction-diffusion equations have been proposed 



IJ. This 



as models for morphogenesis [11311 . pattern formation [g, lZD and self-organization 114 , 
class of equations usually includes the following two features: (i) a nonlinear reaction between 
chemical species describing local production or consumption of the species and (ii) the diffusive 
transport of these species due to density gradients. The simple form of the reaction-diffusion 
equation for a system of A^ species is described by the following set of equations 

dpsix,t) 



dt 



DsApsix, t) + Rs, l< s <N, 



(1) 



where ps{x, t) is the mass density or concentration of species s at time t and location x, A is 
the Laplacian operator with respect to spatial coordinate x, and D^ is the diffusion coefficient of 
individual species s. In this work, we assume that Dg is isotropic and independent of x. The last 
term on the right hand side, Rg, is the reaction term. This term depends on the local density or 
concentration of the individual reacting species and the reaction mechanism governing the system. 
In most pattern forming systems, Rs usually contains non-linear or autocatalytic reaction terms 
with product of the densities of the reacting species. 

Due to their great importance both in biology, environmental science and industry, there has 
been growing interest in a study of these systems both experimentally, by numerical integration 
of the governing equations and via well-tuned analytic theories (see e.g. [|l6l- l22|l and references 
therein). However, solving problems with complex geometry (as is sometimes the case in biolog- 
ical systems) often requires a more efficient and robust method. The Lattice Boltzmann method 
has met significant success in simulating a wide range of phenomena in complex geometries over 
the last decades ll23l - l27n . In contrast with other traditional numerical techniques which only focus 
on the solution of the governing macroscopic equation, the Lattice Boltzmann method is based on 



kinetic theory. In cell-scale modeling of micro-organisms [|28 



30|l for instance, the kinetic nature 



of the Lattice Boltzmann method makes the approach computationally less demanding and allows 
for a relatively simple implementation of microbial interactions between cells. Furthermore, for 



problems involving large domain sizes, the local nature of LB operations allows easier implemen- 
tation on parallel computational platforms thus enabling fast and large scale computations. In 
addition to the above features, the inherent capability of the LB approach in dealing with irregular 



boundaries, makes it suitable for studying reaction-diffusion phenomena in porous media 113111 at 
the pore scale. However the accuracy and efficiency of a numerical method are often evaluated in 
terms of the smallest truncation error within the method. In previous Lattice Boltzmann studies 



of reaction-diffusion equation [1321. I33l1. it is rather unclear as to how the truncation error varies 
with the system parameters such as reaction rate and diffusion constant. These parameters become 
important in pattern forming systems where non-linear reaction terms are present and reaction rate 
as well as diffusion constant may vary over a wide range. Thus, for a better performance and 
accuracy, it is important to find out whether there is a range of optimal parameters that leads to 
the smallest truncation error and a better convergence of the method. Such a study is performed 
in this work for the case of a single species reaction-diffusion systems. Performing a third-order 
Chapman-Enskog multiscale expansion, we investigate the dependence of the truncation error on 
the system parameters. Indeed, for this simple case, while the truncation error linearly varies with 
the reaction rate, it exhibits a pronounced minimum as a function of the diffusion coefficient. 

In order to extend the investigation to a pattern forming multi-species reaction-diffusion model, 
we have selected the Gray-Scott model [34], which serves as a standard paradigm for studying 
reaction-diffusion systems. The Gray-Scott model, though simple, exhibits a wide ran ge o f in- 
teresting dynamical features including spots 1351], spiral waves [36], stationary waves [|37|l and 



spatio-temporal chaos [1380. A particular feature of this model which makes it different from the 



other models is the existence of the so called self replicating spots 1|39D. Spatially localized cell 
like structures grow, deform and make replica of themselves. This act of "cell division" resem- 
bling DNA and RNA replication in cells or the replication growth of biological cells as seen in 
developmental biology makes it an ideal model for studying these biological systems with regard 
to pattern formation. In this reaction-diffusion system, generation of patterns comes usually from 
the instability of an initially uniform state to spatially inhomogeneous perturbations over a certain 
range of wavelengths. The possible range of wavelengths, as determined by a fixed set of system 
parameters, is usually invariant against a change of the system size. A change in system size often 
leads to a corresponding change in the number of spots, stripes or segments observed in the sys- 
tem. Hence, the number of segments or stripes is not invariant but proportional to the system size. 
In contrast, for some biological systems, the pattern forming wavelength is often proportional to 



the system size, while the number of stripes or segments is invariant against the change of system 
size. For instance, some mammalian coat markings have been shown to enlarge in proportion 
to system size ll40n . patterns in some micro organisms like Hydra and Dictyostelium discoideum 



have also been observed to show proportionality with size 114 lh . Modeling this type of biological 



systems with Turing-type reaction-diffusion therefore requires rende ring the governing equations 



dimensionless and adjusting the system parameters in a proper way ll42l l43n. One such approach 



involves using diffusion constants which depend on the concentration of a system size-dependent 



auxiliary chemical factor [|44l-l47h or using the possibility that the concentration of some chemical 



changes with some power of the system size [|48h . Interestingly, it is possible to change the length 



scale of the patterns in the Gray-Scott model via a simple rescaling of all the involved diffusion 
coefficients by the same factor, while keeping all the reaction constants unchanged. We provide a 
test of the validity of this simple approach with Lattice Boltzmann simulations. 

The paper is organized as follows. In the following section, we briefly introduce the Lattice 
Boltzmann simulation scheme for reaction-diffusion equation. We then provide some benchmark 
tests for our LB simulation by comparing our results with analytical solutions for the transfor- 
mation reaction and diffusion of a point source in a domain with periodic boundary conditions. 
Excellent agreement with the analytical solutions is found. We also carry out a truncation error 
analysis of the model via a third-order multiscale expansion. Results obtained from this analy- 
sis are in agreement with our numerical simulations. In section Hill we present the Gray-Scott 
model and, using linear stability analysis, determine the parameter range for the existence of un- 
stable solutions which we identify as a necessary condition for pattern formation. Our numerical 
simulations show good agreement with the predictions obtained from linear stability analysis. In 
section HV] we present a detailed study of the patterns which may be obtained via large amplitude 
perturbations of a linearly stable state. This case comprises the self replicating spots. 

II. THE NUMERICAL MODEL AND ITS VALIDATION 

A. The lattice Boltzmann method 



The Lattice Boltzmann method [|49l- l52|] can be regarded as a mesoscopic particle based nu- 
merical approach allowing to solve fluid-dynamical equations in a certain approximation, which 
(within, e.g. the so called diffusive scaling, i.e. by choosing At = Ax^) becomes exact as the 



grid resolution is progressively increased. The density of the fluid at each lattice site is accounted 
for by a one particle probability distribution fi{x, t), where x is the lattice site, t is the time, and 
the subscript i represents one of the finite velocity vectors ej at each lattice node. The number 
and direction of the velocities are chosen such that the resulting lattice is symmetric so as to easily 
reproduce the isotropy of the fluid [|53ll . During each time step, particles stream along velocity 
vectors Cj to the corresponding neighboring lattice site and collide locally, conserving mass and 
momentum in the process. The LB equation describing propagation and collision of the particles 
is given by: 

fi{x + Ci, t + 1) - fi{x, t) = nji{x, t), (2) 

where f2j is the collision operator. 



The most widely used variant of LB is the lattice BGK model [15 ih . which approximates the 
collision step by a single time relaxation towards a local equilibrium distribution f^'^. The lattice 
BGK model is written as: 

f^'^ix^t) - Mx,t) 



fi{x + ei,t + l) - fi{x,t) 



TLB 



(3) 



/r(^,t)=^,p 



(4) 



where tlb is the relaxation time and the equilibrium distribution //"^ is closely related to the low 
Mach number expansion of the Maxwell velocity distribution given as [|54n 

In Eq. (31), Cs is the sound speed on the lattice (c^ = Ax^/3At^) and Wi is a set of weights 
normalized to unity. The weights Wi depend on the number of velocities used for the lattice. In 
this work we use the two dimensional nine velocities (D2Q9) model with weights Wi given as: 

4/9 e, 
1/9 e, 
1/36 e, 

In order to model the reaction-diffusion equation in the frame work of a Lattice Boltzmann 
BGK model, we introduce a multi-species distribution function /j ^ where the subscript s runs 
over the number of species s = 1...N. In addition, as we are modeling chemical reaction and 
diffusion with no accompanying advection by the solvent velocity, the mean flow velocity u in 
Eq. dH) can be set to zero. This leads to 



Wi 



(0,0),z = 0; 

(±1,0),(0,±1) z = L...4; 
(±l,±l),i = 5....8 



(5) 



/r(^' t) = ujips 



(6) 



Equation ^ satisfies the requirement J2i=o ftl — Ps- The chemical reaction is modeled by in- 
cluding a source term, Rg, in the collision step. This leads to 

fi^s{x + Cj, t + 1) - fi,s{x, t) = h WiRs, (7) 

TlB,s 

where tlb.s is the relaxation time for species s. The source term Rg represent the rate of change 
of density of the species, s, with regard to reaction kinetics. The exact form of the relation between 
the reaction rate Rg and the density (concentration) of each species depends on the type of reaction 
being modeled. The density of the species s, ps, is then computed from the distribution function 

EN r 
i=ofi,s- 

Near equilibrium and in the limit of small Knudsen number (= mean free path/characteristic 
length of problem) the macroscopic reaction-diffusion equation can be recovered using Chapman- 
Enskog multiscale analysis. The detailed analysis leading to the macroscopic equation is outlined 
in the appendix. The relaxation time tlb.s is then found to be related to the diffusion coefficient as 



Ds = c2At(rLB, - 



2r 



B. A transformation reaction 



We provide here a simple test of our Lattice Boltzmann approach for reaction-diffusion systems 
and characterize the truncation error obtained with regard to system parameters. We compute the 
problem of diffusion of a species A undergoing an irreversible transformation or decay reaction to 
a species B, 

A^B. (8) 

The reaction-diffusion equation describing the dynamics of species A can be written as 

— = Da/\Pa{x, y, t) - kbPa[x, y, t), (9) 

where Pa(x, y, t) is the density of species A at point (x, y) and time t. Da is the diffusion coeffi- 
cient of A, Kb is the rate of the transformation reaction and A is the Laplacian operator with respect 
to spatial coordinates x, y. Using the initial condition pAi^o, yo, t = 0) = 6{x — Xo)S{y — yo), 
Fourier transformation of Eq. ^ yields 

^P^^=DAq'pA{q,t)-KBPA{q,t), PA{q,0) = 1, (10) 

at 



where p^(g, t) is the Fourier transform of pAiq, t). Integrating Eq. (flOl) . taking the inverse Fourier 
transform and slightly re-arranging the terms one obtains 

Equation (fTTI) is the analytical solution of the problem posed by Eq. Q on a region infinitely 
extended in space. Setting Ra = —kbPa in the lattice BGK approach introduced in Eq. ^, we 
carried out numerical simulations for comparison with the above analytical solution. We set up 
a two dimensional domain of size L^ = 100 and Ly = 100 lattice units. At time t = 0, we set 
Pa{x = Lx/2,y = Ly/2) = 1 , while p^ = on all other points within the simulation box. 
For the whole region of the domain we initialize the density of species B to zero. We impose a 
periodic boundary condition in x and y directions and set the diffusion coefficient of species A, 
Da = 0.02Ax^/At and that of species B to Db = 0.02Aa;^/At . The resuks obtained are shown 
in Fig. [B Note that at time of order tn ~ LI/Da, the effect of the boundary must be taken into 



account 



[|55|] . Taking this into consideration, we have considered times of order t < L^/Da in 
order to simplify our analysis. A comparison of the density profile of species A obtained at times 
t = 400 and t = 420 from the LB simulation and Eq. (fTTI) is shown in Fig. [TJa) for kb = 0.01. 
Quantitative comparison of the data is done by computing the relative error Ep using the definition 






^ ^ ^.,y.r^,^. .. r^. . .... ^ (12) 



where pA,anix,y) is the density field obtained from the analytical solution in Eq. (fTTI) and 
PA,sim{x,y) is the density field obtained from the simulation. The summation is taken over all 
lattice points in the domain. We obtain a good agreement between the LB simulation and the 
analytical solution in Eq. (fTTI) . The error in this case is less than 1%. We note that increasing 
the value of the reaction rate kb leads to an increase in the relative error as evident by comparing 
Fig.fT](a) for kb = 0.01 and Fig.fT](b) for kb = 0.033. A heuristic explanation for this follows 
from the fact that an increase in kb from 0.01 to 0.033 leads to a 3 times reduction in the time 
scale for reaction t^ = 1/ kb ■ Thus, from a physical point of view, reactions occur three times 
more frequently within the same number of LB iterations. In other words, the relative length 
of the Lattice Boltzmann time step increases with regard to the characteristic reaction time, tR. 
As the LBM is linearly accurate in time, the error arising from the reaction part of LBM should 
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FIG. 1 : Comparison of the analytical solution in Eq. ([TT]) and the Lattice Boltzmann simulation (a) Density 
profiles of species A at t = 400 and t = 420 for kb = 0.01 along the line y = yo (b) The same data as in (a) 
but for Kb = 0.033 (c) Behavior of the relative error in the density of A, Ep, with the diffusion coefficient 
Da at Kb = 0.01 (d) The same data as in (c) but for kb = 0.033. 

exhibit a linear dependence on kb- As will be shown below this behavior is indeed substantiated 
by analytical error estimate and in our numerical simulations. 

On the other hand, in the regime where t/j ^ to, one expects the effect of the reaction rate 
on the error to become less significant. Our numerical simulations indeed confirm this effect. To 
further explore this point, we systematically investigate how the relative error in the concentration, 
Ep, varies with the diffusion coefficient Da and the reaction rate kb- Note that the computation of 
the error for each diffusion coefficient is done at the same physical time corresponding to numer- 
ical times when the spatial extent of diffusion in all directions is the same. This ensures proper 
comparison across the range of different diffusion coefficients. The result of this investigation is 
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FIG. 2: (a) Log-log plot of Ep versus kb for different values of Da- The curves are parallel to the solid 
black line with a slope of 1. (b) Plot of DA,min versus kb- The curve is in line with Eq. (fT4l) . 



shown in Fig.[T](c) and (d) for kb = 0.01 and kb = 0.033 respectively. The first observation from 
both curves in Fig. [T] (c) and (d) is the existence of a minimum in the magnitude of the relative 
error. It is interesting that a similar minimum in the LB truncation error is also observed in the 
computation of shear stress jSq] . Another observation is the shift in the minimum of these curves. 
Figure[I](c) shows that the minimum occurs at DA,min = 0.012, for kb = 0.01, while increasing 
Kb to 0.033 leads to a shift in the minimum to DA,min = 0.03 (see Fig.[T](d)). In order to validate 
the existence of this minimum in this model and characterize the shift observed for an increas- 
ing reaction rate, we perform a third-order Chapman-Enskog expansion of the Lattice Boltzmann 
BGK model for reaction-diffusion equation (see the appendix) and obtain an expression for the 
truncation error up to the third-order in the expansion parameter e. Using Ra = —hbPa in Eq. 
(IA23I) . setting At = 1 and rearranging the terms in orders of ta leads to 



E = 3c%dl^pA ( r^B,A 






(13) 



The third-order LB truncation error, E, is thus the product of time and spatial derivative of density 
with a quadratic polynomial in ta- This polynomial has a minimum at 



Tlb.a- 



^cldtdloA 



+ 0.5 



D 



A,min 



^dtdloA 



(14) 



where we also used the relation between diffusion coefficient and the LB relaxation time. Strictly 
speaking, the value of DA,min not only depends on the reaction rate kb but also on space and time 
variables (through p^). The fact that a minimum does indeed occur in Ep as a function of Da is, 
therefore, not at all a trivial consequence of Eq. (fT3] ). Indeed, the shape of Ep(Da) is not similar 
to a parabola suggesting that the non-trivial effects related to time and spatial derivatives of pA are 
present. 

Nevertheless, it is worth testing to which extent useful information of the behavior of LB trun- 
cation error can be gained via the above analysis. For this purpose, we note two important features, 
which can be extracted from Eqs. [T3]and [141 The first one is that Ep could be a linear function of 
the reaction rate kb (see Eq. (fT3l)). The second observation is that also the position of the minimum 
in Ep{Da), i.e. the value of DA,min could be a linearly increasing function of kb- As illustrated 
in Fig. [21 results of Lattice Boltzmann simulations do confirm the validity of these two aspects. In 
summary, careful choice of the diffusion coefficient and reaction rate can lead to a better accuracy 
of the model. In a multi-species system with many reaction rate constants or wide range of time 
scales, the above discussion may provide guidance in choosing an optimal diffusion coefficient for 
each rate constant. 

III. THE GRAY-SCOTT MODEL 

We consider next the Gray-Scott model as a typical example of a two species reaction-diffusion 
system where the non-linear reaction terms between the species coupled with the transport by 
diffusion give rise to spatio-temporal patterns. The Gray-Scott describes the kinetics of a simple 
autocatalytic reaction in an unstirred homogeneous flow reactor ll34ll . The reactor is confined in a 
narrow space between two porous walls in contact with a reservoir. Substance A whose density is 
kept fixed at Ao in the reservoir outside of the reactor is supplied through the walls into the reactor 
with the volumetric flow rate per unit volume kj. Inside the reactor, A undergoes an autocatalytic 
reaction with an intermediate species B ata rate ki. The species B then undergoes a decay reaction 
to an inert product C at a rate k2. The product C and excess reactants A and B are then removed 
from the reactor at the same flow rate per unit volume kf. The basic reaction steps are summarized 
as follows 

A + 2B^3B, (15) 

B ^C. (16) 



10 



The reaction in Eq. (1151) is the cubic autocatalytic reaction in which two molecules of species B 
produce three molecules of B through interaction with the species A. The presence of B stimulates 
further production of itself, while the presence of A controls the production of B. Substance A 
is sometimes called the inhibitor and B the activator. By constantly feeding the reactor with a 
uniform flow of species A while at the same time removing the product and excess reactants, far 
from equilibrium conditions can be maintained. The equations of chemical kinetics which describe 
the above situations and include the spatio-temporal variations of the concentrations of A and B in 
the reactor take the following form: 

BA 

— = kf {Ao -A)- hB'A + DaV'A, (17) 

— = -[ko + k2)B + k.B^A + DbV^B, (18) 

where A and B are the density of species A and B respectively, Aq is the density of A in the 
reservoir, while Da and Db are the diffusion coefficients of species A and B respectively. In 
order to understand and control the relevant time scale and length scale of patterns observed in the 
system, we introduce variables in the form of time and length scales that represent the physical 
processes acting in the system. For A, the characteristic time scale is the time for the removal of A 
given as 1/kf whereas for B it is l/(/i;/ + /i;2). The characteristic time and length scales for these 
quantities are then: 

TA = l/kf, TB = ll{kf + k2), lA = {DArAf\ lB = {DBrBf^. (19) 

Furthermore, we introduce the following dimensionless quantities: 

A = A/Ao, B = B/Bo, Bo = i-^J . (20) 

Introducing the quantities in Eqs. (fT9l ) and (|2Q|) into Eqs. (fTTI) and (fTSi) leads to: 

BA 
TA^ = -B'A + 1 - i + llV'A, (21) 

CM 

and 

^B^ = +r]B^A -B + llV'B, (22) 

where the parameter ri = —- is the strength of the activation process. It adjusts the 

[kf + k2) 

strength of the non-linear term in Eq. (|22l) . 

11 



The number of parameters can be further reduced by rescaling the time and length scale in units 
of ta and I a respectively. This yields : 

BA 

-^ = -&A + 1-A + V^i, (23) 

dt 



IdB ~^~ ~ 1 ~.~ 

- -- = ^i]B'A -B + -V^B. (24) 

T dt s^ 



where r = ta/tb and s = lA/lB=y'TADA/TBDB- The parameter r describes the relative strength 
of the reaction terms. 

In general, equations (l23l) and (l24l) are difficult to investigate by analytic means. However, 
simple cases exist for which analytical solutions can be found. We start with probably the most 
simple situation of a spatially homogeneous distribution of A and B (V'^A = 0, V^-B = 0). Li 
this case, the steady state solutions of Eqs. (l23l) and (|24)) denoted as A^, and B^ obey 



B^^A, + 1-A, = 0, 



r]BlA, - 5e = 0. (25) 

Equation (l25l) has three solutions. The first solution is the trivial homogeneous solution B^ = 
0,Ae = 1. This state exist for all system parameters. The other two solutions exist provided that 
r] > 2. These are given by: 

At = n^^/tII^ and 5> = !L±v5!Zi. (26) 

A. Test of simulations in the case of spatially homogeneous dynamics 

As a check of our simulation approach, we study the homogenized form of situations where 
Eqs. (l23l) and ((24)) are accessible to an analytical solution. For this purpose, we consider the case 
of spatially homogeneous dynamics with r = 1, implying ta = tb- In this case, multiplying Eq. 
(l23l) by 7] and adding the result to Eq. (|24l) leads to 

d(r]A + B) , ~, f,, dP 

Ai L = ^ _ (^A + 5) ^ -^ = -P (27) 

dt at 
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FIG. 3: (a)A plot of —P{t) = r]A{t) + B{t) — t] versus time for r = ta/tb = 1 but different values 
of the parameters ta and r] as indicated. The solid black curve in each case corresponds to theoretical 
prediction for the parameter values used in the simulation. In all studied cases, the simulation results 
exponentially decay to zero and show a perfect agreement with the analytical prediction, Eq. (1281 ). (b) Plot 
of Logg(P(t)/P(0)) versus I/ta for the same data as shown in (a). 

where P = rjA + B — rj and we used the fact that drj/dt = 0. Equation [27] has the simple solution 
P(t) = P(0)exp(-t) = P(0)exp(-t/rA). In other words, 



r]A{t) + B{t) -r]= [r]A{0) - B{0) - r]] exp(-t/rA) 



(28) 



A test of Eq. (|28l) is provided in Fig. [3] for r = 1 but different values of the parameters ta and 
7]. In the case where r/ < 2, the simulation starts from a spatially homogeneous state A{0) = 1, 
-B(O) = with an additional density fluctuations 5A = 0.5 and 5B = 0.25 added homogeneously 
to A and B, respectively. This is done to break the symmetry which would keep the system at the 
initial state (due to the autocatalytic nature of the Gray-Scott model, without B, no reaction will 
take place). For all other cases where, 77 > 2 we start from the non-trivial states {Af, Bf) given 
by Eq. (|26|) with an additional small fluctuations of the form 5A = 0.1 and SB = 0.1. For all 
values of ta and r] investigated, a perfect agreement is found between theory and simulation. 

B. Stability analysis of spatially homogeneous states 



We proceed in this section to determine the stability of the stationary and homogeneous solu- 
tions obtained in Eq. (|26l) with regard to a spatially homogeneous perturbation. Our analysis starts 
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by looking at the growth rate a of an infinitesimal perturbation about the steady state 

i = ie + 0Ae"*, fi = 4 + 0Be"*, (29) 

where (pA and (ps are the amplitude of the perturbation to the species A and B, respectively. 
Substituting Eq. (|29l) into Eqs. (|23] ) and (l24l) . after linearizing and re- arrangement of the terms one 
arrives at the eigenvalue equation 

(J - aJ) = 0, (30) 

where / is the identity matrix, 4> = i<t>A, 4>bY and the matrix J is given as 

r(2r/i±5± - 1) T^Bf 
J = . (31) 

-2AtBt -{Bf + 1)_ 

The eigenvalue equation in (l30l) has the characteristic polynomial 

Q!^-atrJ+ |J| = 0, (32) 

where tr J and \J\ are the trace and determinant of matrix J . The pair of solutions or eigenvalues 
of matrix J is written 



«i,2 = \ (tr J ± v/(trJ)2-4|J|) . 



(33) 

The eigenvalues ai 2 in Eq. (|33l) . can either be real or complex conjugate depending on the 
relative magnitude and sign of the determinant \J\ and trace trJ. If the real part of at least one 
eigenvalue is positive, the considered solution is unstable. 

For the trivial state {Ae = 1, 5e = 0), trJ = (r + 1) and \J\ = t. Using Eq. ([331) and the 
fact that r > 1 one obtains that both eigenvalues are negative. Hence this state is linearly stable 
with respect to spatially homogeneous perturbations. Next we consider the non-trivial stationary 
homogeneous solutions. In this case, inserting the solutions of A^ and Bf given in Eq. (|26|) in 
Eq. (|3TI) . one obtains that trJ = (t — ^Bf\ and \J\ = t irjBf — 2 j. Furthermore, since in 
this case rj > 2, and using Eq. (|26l ). one can easily verify that \J\{Bf) = {rjBf — 2) < and 
\J\{B~) = {rjB" — 2) > 0. Now from Eq. (l33l) it follows that, independent of the sign of trJ, 
one of the solutions ai 2 is always positive, provided that \J\ < 0. Hence, the state Bf is always 
unstable. For the state B~ , on the other hand, both solutions cti 2 will have the same sign as tr J 
and thus the state B^ may be stable provided trJ < (i.e. r < r]B~). 
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FIG. 4: Plot of the stationary homogeneous state solutions of species A and B given by Eq. (|26] ). Above the 
bifurcation point (r/ = 2.0), two solutions exist: one is unstable to homogeneous perturbations (indicated as 
dashed line) and the other may be stable (plotted as a solid line). At r/ = 2.0, the stable solution switches 
to the trivial homogeneous state (1,0) and for r] < 2.0 only the trivial state exist. The Lattice Boltzmann 
simulation (indicated as symbols ) shows good agreement with the theory. 

Figure m shows a typical bifurcation diagram for the system. We plot the homogeneous steady 
state solutions obtained in Eq. (|26|) as a function of the control parameter 77. For 77 < 2 there exists 
only the trivial state (1,0), while when rj > 2 two additional states {A~, B^) and (A^, B^) emerge. 
The state (A^, B^) is always unstable (indicated by dashed line) while the state {A~ , B~) is stable 
if r < rjB'^ . In the same figure we have also plotted the steady state solutions obtained from the 
Lattice Boltzmann simulation of the spatially homogeneous solutions with small homogeneous 
perturbations at time t = 0. As seen in Fig.lH the LB simulation well reproduces the analytically 
predicted stability diagram. 



C. Inhomogeneous state and Turing instability 

The Gray-Scott model develops a Turing instability for a range of parameters. In this region of 
the parameter space the homogeneous steady state solution becomes unstable and a new stationary 
but inhomogeneous state characterized by the formation of patterns becomes stable. We examine 
the condition for Turing instability in this system by looking at the growth rate a of an infinitesimal 
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spatially inhomogeneous perturbation to the steady state solutions 



A = Ae + 0Ae°*e^^'^ B = 3^ + (pse^'e 



at iqx 



(34) 



As in the case of Eq. (|29l) . 0a and 0b are the amplitude of the perturbations to the species A 
and B respectively, and q is the wave number. Again, inserting Eq. (l34l) into the kinetic Eqs. (l23l) 
and (l24l) and after linearizing and slight re- arrangement one arrives at the eigenvalue equation 



(M -al)4) = 0, 



(35) 



where the matrix M is written as 

r (2riAfB, 



M 



For the trivial solution {A^ 



tBf 



1-^ 

£2 



rr/5f 



-(g^ + i^r + l) 



1, -Be = 0) the matrix M reduces to 





(36) 



M 






(37) 



-ie + 1) 

and solving the eigenvalue equation in (l35l) . we obtain the eigenvalues of the trivial state as 
«! = — r(l + (f /e^) and a2 = — (g^ + 1). Since both eigenvalues are negative, the trivial 
homogeneous state {A^ = 1, i?e = 0) is linearly stable for all system parameters, and independent 
of the wavelength of the applied perturbation. However, it is important to emphasize that this 
stability is restricted to infinitesimal perturbations. Indeed, the trivial state is found to be unstable 
with respect to large amplitude spatially inhomogeneous perturbations. In fact, it is in this regime 
that the so called self replicating spots are observed. The original parameterization of the Gray- 
Scott model by Pearson ll34ll is also based on the numerical simulation of spatially inhomogeneous 
perturbations of the trivial state. 

For the remaining non-trivial solutions, we insert the homogeneous steady state solutions A^ = 
l/(r]Bf) into the matrix M in Eq. (|36l ) and solve for the eigenvalues of M with the characteristic 
Eq. (|32l) by replacing J with M(q). The corresponding eigenvalues are then obtained from 
Eq. (133]) by replacing trJ and \J\ with tiM{q^) = r - r]Bf - q^ir/e^ + 1) and \M{q'^)\ = 
(f'T/e^ + (f{Tr]Bf/e'^ "''")+ T~{lBf — 2) respectively. 
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FIG. 5: The plot of |M"(g^)| versus (f for states (a)i3+ (b)Se • These plots show the range of growth modes 
q for which the determinant |M"(g^)| is negative with the possibility of pattern formation. 

Turing structures or patterns emerge when the system becomes unstable with respect to in- 
homogeneous perturbations. Again, at least one of the eigenvalues becomes positive (unstable), 
when \M{q^) | < 0. \M{q^) \ is a parabola in q^ which attains its minimum value for 






{e' - VBI 



(38) 



Since g^ > 0,aminimumin |iVf(g^)| exists only if e^ > rjEf. This is one of the conditions for 
Turing instability in this system. The boundary of the instability band or range of wave number q 
for which \M{q^) | < is given by the roots of the equation \M{q^) | = : 



9?,2 



{^Bf - e') ± JivBt - e'Y - Ae^^Bf - 2) 



(39) 



In Eq. (1391) . there is an important observation concerning the state 5^ . Using rj > 2 and Eq. (|26l) 
one can show that for the state B^, the condition rjB^ < 2 holds for all ?7 > 2. The consequence 
is that Eq. (|26|) has always one negative root and one positive root independent of the value of e'^. 
Since \M{q'^)\ < for g = (see Fig. I2a)), this opens already an instability band for pattern 
formation as regards the state B^. On the other hand, for the state B~, two distinct positive real 
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FIG. 6: The phase diagram of the model at Da/Db = 2 showing the Turing curve and the Hopf curve for 
parameter space spanned by r and t]. For r] < 2 only the trivial state exists and the self replicating spots 
(SRP) are observed in this regime only. For r/ > 2 two additional states {Af, Bf) emerge. In this figure we 
have only shown the states {A~ ,B~). The majority of the patterns formed from the B~ state in our Lattice 
Boltzmann simulations are observed in the shaded region of the Turing space. 

roots are necessary for an instability band of patterns. Thus, the following condition has to be 
satisfied 

r]B~ < ^2 < -8 + 11.65?75~ (Turing space). (40) 

The first condition, e^ > rjB" as discussed above, is necessary for the formation of Turing patterns 
while the second one reflects the requirement of a positive discriminant in Eq. (|39| ). Using the 
definition of e^ the first condition can be re- written as 



/2 



TF > ^b;. 



(41) 



This means that the diffusive length scale for the species A {I a = ^/Data) must be at least rjBl 



times larger than that of B {Ib = \/Dbtb)- In other words, for a given value of parameter r, 
the diffusion coefficient of species A has to be [rjB" /t) times larger than that of B. In 77 and r 
parameter space, the curves r = {Db/ DA)riB^ , and r = (Db/Da){—8 + 11.65r]B^) define the 
limits of stability with respect to Turing patterns. The first curve r = {DB/DA)f]B^ is plotted as 
the Turing curve in Fig.[6]for Da/Db = 2. The second curve r = [Db/ Da){—'& + ll.QhrjB') 
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FIG. 7: (a) Spatial distribution of the density B at time t = 400000, r? = 2.0139 and r = 2.7330. The for- 
mation of stripes can be observed at these parameters, (b) Amplitude of the Fourier components of density 
fluctuation {B - Bg) at time t = 400000, as function of the wave vector g = 27r (^{ux/Lx)'^ + {uy/Lyf') . 
The dimensionless density Bf. corresponds to the unstable homogeneous state for the selected set of param- 
eters Tj and e. The white region in the Fourier spectrum corresponds to the excited wave numbers {ux, ny). 

lies above the first curve and falls outside the plotted range. It is therefore not shown in the figure. 
At zero mode (g = 0), another important instability known as the Hopf instability occurs when the 
real part of a pair of complex eigenvalues passes through zero. In other words, a Hopf instability 
characterizes the transition from a decaying oscillating mode (triVf (0) < 0) to an oscillation with 
growing amplitude (trAi"(0) > 0). Thus, the limit of Hopf instability is given by the condition 
trA4"(0) = 0. The dashed black line in Fig. Vindicates the limit of the Hopf instability. The small 
dashed area in the Turing space is the region where most patterns are expected to be observed. 

D. Lattice Boltzmann simulation of the spatially inhomogeneous dynamics 

In this section we perform Lattice Boltzmann simulations of the Gray-Scott model for different 
values of parameter r and rj. Here our simulation and parameterization is based on the non- 
trivial states (y4^, Bf). Starting from the homogeneous steady state (A^, B'^), we apply a small 
amplitude density fluctuations of the form 5p = (f)cos{qxx) cos{qyy) where is the amplitude 
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FIG. 8: Stable time independent Turing structures developed from infinitesimal perturbations to state 
{A^,B~) at pai-ameters (a) i] = 2.0072819, r = 2.742424 (b) r? = 2.013958, r = 2.73333 (c) 
r] = 2.017971, r = 2.707462 and (d) ?/ = 2.0184336, r = 2.7272. The system size in all cases con- 
sidered above is 200 x 200 lattice units. 

and q is the wave number of the perturbations. We have chosen = 0.001, and Qx = Qy = ^ 
in the simulations. Fig. ITfa) shows a developed stable structure from the small amplitude initial 
perturbation to the B^ state with parameters r] = 2.014 and r = 2.733. For the purpose of 
comparison with the prediction of linear stability analysis, we perform the Fourier transform of 
the pattern in Fig. 13 a) and calculate the excited wave numbers in the Fourier spectrum using the 
relation 

q = 27r {{njLxy + (n./L,)^)'/' , (42) 

where n.^ and Uy satisfy —Lx/2 < n^ < L.J./2 and —Ly/2 < Uy < Ly/2 respectively. 
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The Fourier spectrum is shown in gray scale in Fig. |7]^b) and the excited wave numbers are 
{nx,ny) E {(±1, ±5); (±3, ±4); (±4, ±3); (±5, ±1)}. Using these values of n^ and riy in Eq. 
(I42l) to calculate q and the parameters t] = 2.0139 and e"^ = 5.466 in Eq. (|39l) . we found that 
all the excited wave numbers from the simulation fall within the instability band predicted by 
linear stability analysis in Eq. (|39l ). This provides a further validation of our Lattice Boltzmann 
simulation with regard to this model. 

By performing a number of similar simulations with different values of r and t] we have found 
that Turing patterns develop over some part of the region where the B~ state is Hopf or Turing 
unstable (indicated as the Turing space in Fig.|6l). The panels in Fig. [8] show developed stationary 
structures for typical values of the parameters t] and r. These parameter values fall between the 
saddle node bifurcation curve (rj = 2) and the Turing curve (see Fig. (6]). One observation from Fig. 
[8]is that increasing the value of t] within the Turing regime leads to the development of a lace-like 
structure in the patterns. 

IV. BEYOND LINEAR STABILITY: SELF REPLICATING SPOTS 

In this section, a further example is provided for the maturity of the Lattice Boltzmann method 
in studying pattern formation within the Gray-Scott model. The patterns discussed so far are in 
the part of the phase diagram (Fig. |6l), where linear stability analysis predicts that homogeneous 
solutions are unstable with respect to small perturbations. There are also other types of structures, 
which occur in a regime, where the trivial homogeneous state is linearly stable. These patterns 
emerge only if the homogeneous state is perturbed strongly enough. A prominent example of this 
type of structures are the so-called self replicating spots. 

Figure|9]illustrates the patterns emerging from a finite amplitude perturbation of the trivial state. 
The starting configuration corresponds to a rectangular box of species A and B with densities 
A = 0.5, B = 0.25 placed at the center of a domain filled with species A and B at densities 
A = LO, B = 0, respectively. Obviously, such an initial state represents a strong perturbation 
of the trivial state (A = 1, B = 0, i.e. A = Aq, B = 0). The sequence of images in Fig. |9] 
demonstrates how spots form, elongate and then replicate as time proceeds. This self replication 
process continues until the whole simulation cell is filled with the spots. Interestingly, the number 
of spots increases with time, while the size of an individual spot seems to remain roughly constant. 
We have repeated this simulation for a larger box size but otherwise exactly the same parameters. 
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The result of this study is also shown in Fig. [9l As seen from the last image in Fig. |9l the size of a 
spot does not change, but only the number of spots increases as to fill the entire simulation cell. 

In the light of above presented results, one may rise the question, whether it is possible to keep 
the number of spots constant but tune their size. An answer to this question is obtained by noting 
that the dynamics of Eqs. ( |23l) and (l24l) depends only on the dimensionless quantities e, r and r]. In 
other words, we must check, whether it is possible to tune the length scale of the problem without 
altering the values of these dimensionless parameters. This would ensure that the thus obtained 
new solution will have exactly the same shape (and thus the same number of spots) but a different 
length scale (different size of spots). Indeed, a look at the parameter e reveals that it is equal to 



the ratio of two characteristic lengths Ia and Ib, s = lA/lB=yTADA/TBDB- Thus, if we multiply 
both I A and Ib by a constant factor A, the parameter e remains unchanged. Furthermore, in order to 
keep also the other two parameters r and t] constant, the simplest choice to achieve such a change 
of length scale is via diffusion coefficient, i.e. via Da -^ ^^Da and Db —^ ^^Db- 

In order to test the above idea, we design two systems such that the system 1 has a linear 
dimension of Li = 200 lattice units with diffusion coefficients Da,i = 0.016Ax^/At and Db,i = 
0.008Ax^/At. For the system 2, we choose L2 = 400 lattice units, which means that A = 2. 
Following the above arguments, we set the diffusion coefficients of the species A and B in the 
system 2 to Da,2 = A^/^^.i = 0.064AxVAt, and Db,2 = \'^Db,i = 0.032AxVAt, respectively. 
As initial state, we perturb the trivial state exactly in the same way as described in the context of 
Fig. m and impose periodic boundary conditions in both the x and y directions. Note that the size 
of the square perturbation must also be multiplied by A in conformity with the change of length 
scale. Results of these simulations are shown in Fig. [TOla) and (c). The structure of the patterns 
is identical for the two systems within numerical discretization errors. The inner core diameter of 
the spots is found to scale as the diffusion length of species A, I a- A more quantitative comparison 
of the data is provided in Fig. [TOlb) and (d), where time evolution of the density profiles is shown 
for the both studied system sizes in a space-time plot along the x direction. 

The above arguments on how to tune the length scale while keeping the shape of the patterns 
unchanged is quite general and applies to any other solution of the Gray-Scott model as well. Here, 
we provide an example from the Turing regime. This is an interesting test, as linear stability anal- 
ysis predicts that when the system size is increased, the number of stripes or segments is increased 
accordingly. However, this applies only if all other parameters are kept constant. Interestingly, by 
proper regulation of the diffusion coefficient, our numerical simulations in the Turing regime con- 
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FIG. 9: Snapshots of the density distribution of species B showing different stages in the self replication 
process at times (a) t = (b) t = 10000 (c) t = 20000 (d) t = 50000 (e) t = 100000 (f) t = 300000 with size 
200 X 200 lattice units and (g) t = 300000 with size 400 x 400 lattice units. 
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FIG. 10: Snapshot of the spatial distribution of the density B showing self replicating spots at time 300000, 
r] = 1.86, T = 3.40 and e = 2.61 for the two systems with lattice size (a) 200 x 200 lattice units (c) 400 x 
400 lattice units, (b) and (d) show space time plots of the density profile of the self replicating spots along 
a line in the y direction for the two system size in (a) and (c) respectively. 

firm that it is possible to make the wavelength proportional to the system size and keep the number 
of stripes or segments invariant. Results of these simulations are shown in Fig.fTT] For parameters 
r] and r in the Turing regime, we choose r] = 2.014 and r = 2.733 and consider two systems with 
a scaling factor A = 2. The spatial density distribution obtained from the simulations is shown 
in Fig. [rTT a) and (b) for the two systems respectively. Not unexpectedly, the patterns exhibit the 
same structure with equal number of stripes and segments. To further support this observation, we 
carried out numerical simulations over a range of system sizes from 50 to 500 lattice units. The 
thus obtained density profiles along the line y/L = 0.85 are plotted in Fig. \T2\a.) for all studied 
system sizes. For the sake of visibility, each individual curve is shifted by a multiple of 2 along the 
vertical axis. It is clear from the figure that the number of stripes does not change with the system 
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FIG. 11: Turing pattern showing spatial distribution of the density B~ at time t = 400000, rj = 2.016933 
and r = 2.73030 for two systems of size (a)200 x 200 lattice units (b) 400 x 400 lattice units. The patterns 
are cleaiiy identical in the two cases. 

size. To further emphasize the similarity of the patterns, we directly compare on the same figure 
all the data using the same shift for all the curves. Clearly, the data collapse into a single curve. 

As additional demonstration of wavelength regulation and proportion preservation of the pat- 
terns, we perform Fourier transform of the patterns obtained from each system size in the range 
of 50 to 500 lattice units. We calculate the maximum excited wave number q^ax in the Fourier 
spectra of the density field using Eq. (|42|) . Fig. [T2lb) shows the plot of the wavenumber excited 
with the system size for different diffusion coefficients. It is clear from the plot that the maximum 
excited wavenumber qmax decreases in proportion to L and thus the generated pattern is expected 
to preserve the proportion as observed in our numerical simulation. 

V. SUMMARY 

In this work, we study reaction-diffusion systems via Lattice Boltzmann computer simulations. 
Starting from the analytical solution of a simple prototypical model (a single species undergoing 
transformation reaction and diffusion), we perform a systematic study of the Lattice Boltzmann 
truncation error of the model. We uncover interesting behavior of the truncation error with the 
system parameters. The error is found to have a minimum at a given value of diffusion coefficient. 
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FIG. 12: (a) Density profile of B~ along the line yjL = 0.85, obtained from the Turing pattern in Fig. 
[TTlat time t = 400000, with parameters t] = 2.014 and r = 2.733 pai^ameters . The number of stripes are 
invariant over an appreciable range of system size, (b) Characteristic wave number of the Turing pattern in 
Fig. [TT] plotted against the system size for different values of e^ (realized via a variation of Da/Db)- The 
curve shows preservation of proportionality between the wave number and the system size. 

The position of minimum is shifted for increasing values of the reaction rate constant. These 
observations are in agreement with the analytical findings from a third-order Chapman-Enskog 
multiscale expansion. 

A study of the Gray-Scott reaction-diffusion model is also provided. Here, we perform a linear 
stability analysis of the model and determine the relevant parameter range for pattern formation. 
Lattice Boltzmann simulations of this interesting reaction-diffusion system are found to be in good 
agreement with the predictions of the linear stability analysis. In addition to a test of the linear 
stability phase diagram. Lattice Boltzmann simulations provide valuable information on the details 
of the patterns formed in different regions of the parameter space. An example is the formation 
of striped patterns in most parts of the Turing regime above the Hopf bifurcation curve. Another 
very interesting example is provided by the so called self replicating spots, which lie beyond the 
linear stability regime. Self replicating spots occur via large amplitude perturbations of the trivial 
homogeneous solution. 

Furthermore, a survey of the parameters entering the scale invariant form of the Gray-Scott 
model suggests that the simplest choice to tune the length scale of the obtained patterns (while 
keeping its shape unchanged) is to multiply all the relevant diffusion coefficients with the same 



26 



constant factor without any modification of the reaction rates. Interestingly, this lets the time scale 
of the process unaffected. In other words, in systems with different diffusion coefficients but same 
reaction rates, patterns exhibit the same shape (but different sizes) exactly for the same physical 
time. Results obtained via Lattice Boltzmann simulations confirm this behavior. It is noteworthy 
that this act of regulating diffusion coefficient of species or morphogens as the case may be, is also 
observed in some biological systems. This observation is by no means limited to this model, the 
analysis can also be extended to other reaction-diffusion model. 
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Appendix A: Chapman-Enskog procedure for reaction-diffusion equation 

In this appendix, we derive the macroscopic reaction-diffusion equation from the Lattice Boltz- 
mann model. The LB equation for reaction-diffusion equation is written as 

fUx + e,At, At + t) - fUx, t) = ' ' ^ + MwA (Al) 

TlB.s 

To obtain a corresponding macroscopic partial differential equation from the finite difference Eq. 
(lAll) . we perform a Taylor series expansion of the left hand side of Eq. (lAll) and obtain 

1 n\ tlb.s 

71=1 ' 

The Chapman Enskog procedure introduces two time scales, a fast time scale, ti, associated with 
convective transport and a slow time scale, t2, associated with diffusion. The time derivative is 
then expanded as 

'H = tut -r t L/f , 



dt = edi'^ + e^di^\ (A3) 



The spatial derivative is written as 



d.^=edi'^. (A4) 
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The equilibrium distribution as reaction term are expanded as 



R, = Rf + ei?(i) + e^Rf^ + e'Rf + 0{e') (A6) 



Inserting Eqs. (|A6l ). (|A5] ). (lA4l) and ( |A3] ) in Eq. (lA2l) . one obtains 



\t{ed^^^ + e'df^ + ee.„C) + ^(^^^i'^^f ^ + 2e^e.„af ^ai^) + 6^e.„e.,C^S + 2e'e..9l^)C 



= - (l7'si^,t) - (4? + 6/« + eVf; + eVff + 0{e')))+Atw.{eR^^+e'Rf^+e'Rf^+Oie% 

(A7) 
Grouping terms of the same order in e yields the following successive approximations. 

0(e(°)) : ff^ = f:;i, implying that R<f>^ = 0. 

(A8) 
Note that this condition follows directly from the conservation of mass. 

0{e') : At (a« + e..C) 4? = -^f^ + ^t^^^^' (^9) 

On the level of e, there is no mass diffusion, diffusion process takes place on the scale of e^. Fur- 
thermore, for diffusion driven reactions, the diffusive flux must bring the species together before 
reaction and reaction becomes a second order effect. Since on the scale of e, there is no mass 
diffusion, that implies Rs =0. We would consider this case in this derivation and Eq. (IA9I) 
becomes 

0{e') : At(ame..C)4? = "^/ff- (AlO) 

^ (d^ + 2e..9«C + e.«e.,C5S) 4? = '^f' + ^tw.R^' ■ (AH) 
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O(e^) : At {d?^f!S + d?^fi:^ + (a;^) + e.„C) 4?) 



^(dr+2e..di^'d2^+e..e.,d^d2j)fQ 



(A12) 
Putting the expression for f- J from Eq. (lAlOl) into Eq. (lAllI) yields 

^/f ^ = -At9i^)4? + At^ (r.B. - (5f ) + e..C)' 4? + At^^./^f • (A13) 

In Eq. (IA12I) . we insert the expression for f^J and f^J from Eqs. (lAlOl) and (lAllI) and obtain 

-4? = -A^^f ^4? + At^(2n^B. - 1) (ai^) + e.„C) a;^^4? 



Tlb„ 



-At^ r2 



LB. - -LB. + Ij (a« + e.„aW)' 4?-rLB.At^ (9f ) + e^^aW) ^.i?f )+At^.i?(^). 

(A14) 

Next we take the moments of the distribution functions in Eqs. (IA9I ). (IA13I) and (IA14I) . Note 
that in order to preserve the isotropy of the lattice tensors, the chosen lattice speeds and weights in 
the equilibrium distribution function must obey the following moments or symmetry conditions. 

{a)'^Wi = 1 

i 

{b)'^Wieia = 

i 

(c) ^ WjCiaei/j = cl6ap 

i 

{d)^^Wieiaeipei^ = 

i 

{e)^Wieiaeipei^eis = c'^si^o.^S^s + ^a^Sps + Sas^^y) 

(A15) 

Using Eq. (IA15I) and given that the local equilibrium takes the form f^'^ = fl = WiPs, we impose 
the following conditions of conservation of mass on the equilibrium distribution function 

Yl 4? = p^^ Yl ^^"4? = 0' Yl ^io^^ii^fS = Ps(^%^ (A16) 

i i i 
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We further assume that higher order corrections of the equilibrium distribution do not contribute 
to the local values of the mass, whereby obtaining 

^/\(") = 0, for n>l. (A17) 

i 

Taking Y.i of Eq. dMl) and using Eqs. (IaTSI) . (IaT61) and (IaTtI) yields 

5j'V. = 0, (A18) 

Taking Y.i Eq- (lAlSl) and using Eqs. (IaBI) . (IaT61) . (IaTtI) and (IaTSI) leads to 

Taking ^. of Eq. (lAT4l) and using Eqs. (IATSI) . (IAT61) . (IaTtI) . and (IaTSI) leads to 

9f Vs = -3c^At^ (r,^,,, - TLB, + ^) af )C4>.5„, - r^^.^ed^Rf^ + i?f ). (A20) 

We multiply Eq. (IA18I) by e, Eq. ( IA19I) by e^ and Eq. (IA20I) by e^ and add all this together, thus 
arriving at 

dtPs = clAt Uji,s - ^ j dl^ps + Rs- SAt^cl (tI^^ - tlb.s + ^ j dtdl^ps - Tui^s^t'dtRs. 

(A21) 
We can further re-write Eq. (|A21|) as the macroscopic reaction-diffusion equation and a third- 
order truncation error term E 

dtPs = Dsdl^p, + Rs-E {All) 

where the diffusion coefficient is given by Ds = cj.At (tlb.s — 0.5) and the error term takes the 
form 

E = 3At^cl Tr^B,. - ^LB,s + I) dtdlPs + r^B^edtR, + ^(e^). (A23) 
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